!
! Résolution de Au = g par la méthode de Gauss-Seidel
!
SUBROUTINE GAUSSSEIDEL(n,p,eps,h,alpha,beta,c,seuilResidu,seuilCompteur,u)
!
implicit none
!
! Déclaration des variables d'entrée
DOUBLE PRECISION, INTENT(IN) :: eps, h, alpha, beta, c
DOUBLE PRECISION, DIMENSION(n*n) :: u
INTEGER, INTENT(IN) :: n,p
!
! Déclaration des variables locales
DOUBLE PRECISION :: somme, seuilResidu, residu
INTEGER :: compteur, seuilCompteur, i
DOUBLE PRECISION, DIMENSION(n*n) :: g, vect_residu
DOUBLE PRECISION, DIMENSION(5) :: A
! Duree de programmation
REAL :: t_debut, t_fin
! On lance le chronometre
CALL CPU_TIME(t_debut)
!
! Appel des differentes subroutines
CALL createG(n,eps,alpha,beta,c,h,g)
CALL createA(eps,alpha,beta,c,h,A)
print *, ''
!
! Premier calcul du résidu
u = 0
CALL produitmatricevecteur(A,u,n,p,vect_residu)
vect_residu = g - vect_residu
CALL PRODUITSCALAIRE(vect_residu,vect_residu,n,residu)
residu = DSQRT(residu)
!
compteur = 0
DO WHILE (residu .GT. seuilResidu .AND. compteur .LT. seuilCompteur)
	!
	! Calcul de u
	DO i = 1,n*n
		IF (i.LE.n) THEN
			somme = A(4)*u(i-1) + A(2)*u(i+1) + A(1)*u(i+n)
		ELSE IF (i.GE.(n*n-n)) THEN
			somme = A(5)*u(i-n) + A(4)*u(i-1) + A(2)*u(i+1)
		ELSE
			somme = A(5)*u(i-n) + A(4)*u(i-1) + A(2)*u(i+1) + A(1)*u(i+n)
		END IF
		u(i) = (g(i) - somme)/A(3)
	END DO
	!
	! Calcul du résidu
	CALL produitmatricevecteur(A,u,n,p,vect_residu)
	vect_residu = g - vect_residu
	CALL PRODUITSCALAIRE(vect_residu,vect_residu,n,residu)
	residu = DSQRT(residu)
	!
	compteur = compteur + 1
END DO
! Arret du chronometre
CALL CPU_TIME(t_fin)
!
! Affichage du code de sortie
IF (compteur .LT. seuilCompteur) THEN
	PRINT *, 'Code de sortie : convergence'
ELSE
	PRINT *, 'Code de sortie : seuil de compteur depasse'
END IF
!
! Affichage des résultats
WRITE(*,*) 'Nombre d iterations pour la methode Gauss-Seidel :', compteur
WRITE(*,*) 'Valeur du residu :', residu
WRITE(*,*) 'Temps d execution :', t_fin - t_debut, 'secondes'
!
END SUBROUTINE
